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Electrostatic gyrokinetic absolute equilibria with continuum velocity field are obtained 
through the partition function and through the Green function of the functional integral. 
The new results justify and explain the prescription for quantization/discretization or taking 
the continuum limit of velocity. The mistakes in the Appendix D of our earlier work 
[J.-Z. Zhu and G. W. Hammett, Phys. Plasmas 17, 122307 (2010)] are explained and 
corrected. If the lattice spacing for discretizing velocity is big enough, all the invariants 
could concentrate at the lowest Fourier modes in a negative-temperature state, which might 
indicate a possible variation of the dual cascade picture in 2D plasma turbulence. 
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I. INTRODUCTION 



The statistical equilibrium states of a Galerkin truncated system are called the absolute equi- 
libria. Such equilibria have been used to study the statistical dynamics of various systems since 
Lee, 1 including models for the Bose-Einstein condensates and superflows (see Ref. 2 and refer- 
ences therein,) besides other systems such as classical (magneto)hydrodynamics, drift waves and 
(gyro)kinetic fluctuations in plasmas (see Ref. 3 and references therein), among others. The novel 
feature in (gyro)kinetics is the appearance of the extra velocity variable in distribution functions. 
Different treatments of the extra dimension of the velocity variable would lead to different abso- 
lute equilibria. This paper studies gyrokinetics statistical equilibria in the case without additional 
treatment of the velocity variable. This system is thus closer to the original one. Our focus is on the 
general feature of the problem and the specific solving techniques, with discussions of the results 
in the common context of turbulence physics. As in the article by Krommes, 4 the results presented 
here are of potential interest also to non-plasma physicists. (A minimization of references is thus 
managed to keep the general readers in an easy environment.) This work is self-contained and, in 
particular, knowledge of the results in the paper by Zhu and Hammett 3 (henceforth referred to as 
I) is not required. As supplementary material to further highlight discretization schemes, we will 
actually provide an alternative approach, the partition function calculation, whose intermediate 
steps reproduce the results in I. 

Zhu and Hammett 3 obtained gyrokinetic absolute equilibria which can be considered as an 
approximation of the "exact" continuum result. It is fundamental to ask what the continuum limits 
or the exact results are and how the continuum limits should be taken (this latter issue should have 
been clear from the discretization scheme, but there are subtleties). A deeper understanding of 
gyrokinetic absolute equilibrium is promising for plasma turbulence, as indicated by theoretical 



work in hydrodynamics. For example, Ref. 5 found that in some cases two dimensional (2D) 
turbulence is actually not far from absolute equilibrium, and dissipative three dimensional (3D) 
turbulence is also partially thermalized (see, e.g., Refs. Ill andla) 

As we will see, the problem itself is formulated in terms of integrals over the velocity field. 
Therefore, alternatively, we can perform the calculations without starting with the discretization 
of velocity used in I, and the results should be the continuum limits I. This is the main technical 
contribution of this paper, as detailed in Sec. [TTTJ In other words, instead of first applying discretiza- 
tion and then taking the continuum limit, we will take the strategy of starting from the calculation 



with continuum velocity and then discretize the results. Such an approach provides an alternative 
view of the problem and may clarify the issue better: Suppose we have somehow arrived at the 
correct continuum limit, then figuring out the continuum limit directly will of course justify and 
explain the discretization/quantization scheme and the prescription of taking the continuum limit 
from the discretized results. Also, calculations with arbitrary velocity field offer direct compar- 
isons and allow us to elucidate other issues such as the physical effects of numerical errors and of 
finite truncation of the velocity scale (typically from the collision operator). As explicitly stated 
in the Footnotes 46 and 65 of I, the coauthors disagree on what the continuum limit should be and 
on how to obtain the limit in the 2D case. To avoid confusion and extra burden of references to 
the general audience, we make this paper sufficiently self-contained. We clearly, though briefly, 
point out the mistakes in Appendix D of I. Together with the discretized- velocity results, the con- 
tinuum limits may also promote more fundamental considerations on the statistical dynamics of 
magnetized plasmas. 

H. THE MODEL 

Gyrokinetics is a powerful model for modern research on microturbulence in magnetized plas- 
mas with many applications to both astrophysics and fusion (see, e.g., Krommes 4 - which is in- 
tended also for non-plasma physicists.) Here, we limit ourselves to the nonlinear electrostatic 
gyrokinetic Vlasov-Poisson equations as in I, where more literature can be found. 

The gyrokinetic equation for electrostatic fluctuations in slab geometry with uniform back- 
ground magnetic field B = B z reads 



where we follow the notations and normalization conventions used earlier in I: g(R, vn , v±) defined 
at the gyrocenter is a component of the fluctuating distribution /, where / = exp{— q(p}F + g + 
(v 9 )r + h.o.t., which gives the deviation from the maxwellian F . Here h.o.t. denotes the higher- 
order terms in the gyrokinetic ordering and (-)r the gyroaverage around R, i.e., the average of any 
quantity \I/ along a ring of gyroradius p surrounding the guiding center R and perpendicular (_L) 
to the magnetic field direction (||): 




(1) 



Jv&(r^(rn-R||)^|r ± -R ± |-p(R)]^r 
J6(r\\ -R||)5[|r ± -RJ - p(R)]d 3 r 
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with qp x B = mvi, which illustrates how p and v_l experience the gyroangle together. [Through- 
out this article, unless explicitly specified, J is used to denote the definite integral over the entire 
domain. Using a Fourier representation \I/(r) = J^ k exp(— ik • r)^ k , and considering a straight 
magnetic field for the sake of simplicity, we have (^)r = exp(— ik ■ R)J (kj_p)^/ k , where J 
is a Bessel function. 7 ] The gyrokinetic equation expresses how guiding centers evolve in time due 
to parallel motion along the magnetic field, the gyro-averaged E x B drift across the magnetic field 
(this is the nonlinear term), and the parallel electric field acceleration. This equation is closed by 
using the gyrokinetic quasi-neutrality equation to determine the electrostatic potential, which in 
Fourier space in normalized units is given by 

B(k) f POO POD 

0{k,t) = / d 3 vJ (k ± v 1 _)g = B{k) / dv\\ / dv±v±J {k ± v ± )g. 

J J — oo J 



The physical (dimensional) variables having subscript 'p' are normalized as follows: t = t p v±/L 9 

T V 3 L 

x = Xp/pth, y = Vp/pth, z = z p /L,v = v p /v th , = y? P %> h = and F = F 0p vl/n . The 

equilibrium density and temperature of the species of interest are n and T ; the thermal velocity 



is t> t h = A/T /m; the Larmor radius is p th = v&/Q c where the Larmor (cyclotron) frequency 
is f2 c = qB/m. L is the reference scale length (i.e., the system size), satisfying p/L 1 for 
consistency with the gyrokinetic ordering. 

For plasmas in a two dimensional (2D) cyclic box, with vn being integrated out and the subscript 
x in v± omitted, the whole system in wavenumber space is 

d t g(k, v)-zx P^o(H<£(p) ■ q^(q. v) = (2) 

p+q=k 

and 

0(k)=/3(k) [ vdvJ (kv)g(k,v). (3) 



Here B(k) = r+1 2?r p ^ and f (x) = I (x 2 )e~ x is an exponentially-scaled modified Bessel function, 
where Iq(x) = Jo(ix) and r represents the shielding by the other species which is treated as having 
a Boltzmann response of some form (an anisotropic response model, such as those respecting the 
zonal modes, can also be considered as in I). 

We should remark that the electrostatic gyrokinetic equations constitute an integro-differential 
system. When they are discretized, we denote the corresponding solution by g to distinguish it 
from the exact solution g of the original equations. In general, g(v ) ^ g(v). As our interests are in 
the statistical properties of the Fourier Galerkin truncated system for which g would be reserved in 



the case with continuum velocity, g is used for the solution to the system with further discretization 
of the velocity. To our best knowledge, there is no rigorous mathematical result on the singularities 
yet, not to mention the collisionless dissipation of the gyrokinetic equations. However, we believe 
the Fourier Galerkin truncation will also smooth the velocity space structure and g here should not 
produce dissipation, unlike the shock solutions in Burgers equation, so that it should make sense 
to talk about a statistical equilibrium state. 



III. CALCULATIONS 

The Fourier Galerkin truncation here is defined by setting all Fourier modes beyond the wave 
number set K = {k : k min < \k\ < k max ] (the summation over which will be denoted by Yl) 
to zero. As Lee 1 did for Euler equation for ideal fluid flow, we first observe that the dynamics of 
the "gas", composed of the real and imaginary parts of the Fourier modes (denoted by a), satisfies 
Liouville theorem. Actually 1^ = 0, where the dot represents time derivative. Then, we proceed to 
find the canonical distribution. Let us start with the 2D case. The only known rugged (conserved 
after Fourier Galerkin truncation) invariants for the 2D gyrokinetic system are 

G(v) = /|£/ = £^(M)| 2 
and the mean effective electrostatic potential "energy" 

*=/;>+^ 2 -^H£^(k)i*, 

(c.f., e.g., I and references therein) where V is the volume (area) of the integration domain. Notice 
that G(v) is a function of v, so we have now "one plus a continuum" of conserved quantities. 



A. Green function 

The Gibbs distribution 8 reads 

Z _1 exp{— S}, Z— jDaexp{—S}, 

with 

S = J a(v )G(v )dv + «oE 

and with a running over the configuration space spanned by the real and imaginary parts of g,- 
which we don't expose explicitly. Here the constant of motion S is formed by introducing a and 



a(v), the "(inverse) temperature parameters" as the Lagrange multipliers when maximizing the 
Gibbs ensemble entropy: 

S = t { I n (V ) | </(k. r)\~d c -\- < 1 1)2 7T - ) ( A' ) / r(IrJ it (kr)<i(k. r) j r<IrJ i} (kr)ij l (k- '') 




£ / / ^(P,«)M(p,k,u,u)^(k,i;)dudi;, (4) 



' P ^k2 

where "*" denotes complex conjugate. Clearly, M(p, k, u, v) = 5 Pt kC x (k., u, v), with 

C z (k, u, i>) = a(v)5(u — v) + a 2TT /3(k)vuJo(kv ) Jo(ku). (5) 

[We should clarify that G(v ) is introduced as a density over v in the first line of Eq. © and is then 
changed to a density in the u — v plane in the second line with the help of the Dirac delta function: 
Since (|g(k, t>)| 2 ) = (<?(k, t>)g*(k, v)), the variable v can be regarded as appearing once or twice; 
so, a density of v, or, of u and v with u = v (denoted as "IV" and "2V" respectively) could be 
represented depending on the context, which should be kept in mind in case of misconceptions; 
similarly is the discretized case. ("D" is already used to denote the dimension in configuration 
space, so we adopt "V" here for the dimension in velocity space.) The IV and 2V densities, 
though written symbolically the same, acquire different values.] 
With the definition of functional inverse, 

C x (u, x)C(x, v)dx = 5(u— v) = J C(u,x)C I (x,v)dx, (6) 

we find (for details, see a simple approach with careful observations in Appendix I A H and a general 
functional inverse formula which leads to the same result as established in Appendix I A 2|) 

n\ x^-^) 2na (3(k)uJo(uk)vJ (vk) 

O (K,-U,fJ— — — p , — ~ • (7 ) 



a ( v ) a (uMv)[l+2na mr-^dx 
Here 



C(k,u,v) = (g(k,u)g*(k,v)), (8) 

where (■) denotes the statistical average over the Gibbs distribution J Da(-) exp{— S}/Z [A factor 
2 has been eliminated to respect the realizability condition g*(k, u) = g(— k, u)] , is a 2V density 
in the u-v plane. Therefore, given the Dirac delta function in the right hand side of Eq. ©, the 
IV distribution over v of the spectra density (over k) of g 2 (R, v) can be obtained as detailed in 
Appendix |Bj 

G(k,v) = l/a(v). (9) 



This spectrum is independent of k ("equipartition") since the second term of Eq. © does not 
contribute to the line u = v: see Eqs. (IB II) and (|B2I) — the global picture is that G(v), extended 
as 2V later to facilitate the calculation, was introduced in Eq. © as a IV density to which Eq. 
© is reverted. When the IV density is extended to the u-v plane, a Dirac delta function must be 
introduced as in Eqs. © and ©, and, of course the 2V density is then singular for u = v, unlike 
the IV density. The IV result in Eq. © is the same as the case a = 0, because the diagonal 
(u = v) contribution of the a term is zero (as the 2D Lebesgue measure of a line is zero) and 
the nondiagonal (u ^ v) contributions cancel among themselves. We then calculate the spectral 
density D(k) = ^y(|<£(k)| 2 ) of the "energy" E by 



D(k) = Trp(k) / / C(k,u,v)vJ Q (kv)uJ Q (ku)dudv = — y i 2 g- , . (10) 

J J l + a 27Tf3(k)r-!^ 1 dv 

When discretizing, assuming rectangular lattice for simplicity, with particularly 

Si,j I 'rrii 8{u - v) (11) 
we have the covariance density matrix elements 

«m = m, »,» = s f - Z*T%£ w V-> <12) 

a» 1 + a 2ir(3(k) 2^ wfa t 
which is a 2V discrete density. Here Wi(k) = miViJ (kvi), and rrii is the weight of velocity grid 
point Vi. Actually, now rrii = and, 

ai = a(vi)mi, for z = 1, ...,N, (13) 

where V is the total number of grid points. The tilde above a is similar as for g. (Such discretiza- 
tion corresponds to schemes used by current Eulerian gyrokinetic codes and our thermodynamic 
limit is exactly taken with the lattice spacing going to zero.) In I, the same covariance density 
matrix as in Eq. (fT2l) was obtained by using the discretization of Eqs. © and © as the starting 
point. Especially, to facilitate the calculation, Gi was multiplied by 5ij, which, according to the 
current Eq. (fTTI) . transformed the IV density to 2V density: Gi/rrii =^ Gi, with the notations used 
there. Note that now cy, as well as aj, depend on the size Aw, x Avj of the 2V lattice into which 
the distribution curdling on the IV line u = v is dissolved. If a physical example is still required, 
we would mention that Mandelbrot 9 proposed that the turbulent dissipation of a fluid is curdling 
into fractals, which means a Dirac delta function supported by the fractal object (carrier) imbedded 
in the 3D space (or 4D space time.) The introduction of a lattice for discretization facilitates the 
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computation of box counting dimension whose value is actually 1 for our case here (we have a full 
line, not really fractal.) Especially, we have the spectral density of the (i, i)th "cell" 

\ w. w 1 l\ a 2vr/3(A;)w 4 2 a i " 1 1 

which is a 2V density with the two velocity indexes coinciding. As mentioned, the continuum 
limit of c i: j contains a Dirac delta function which is singular for u = v, so this 2V density c^, is 
singular in such a limit. The corresponding IV density can be obtained by following the procedure 
given in Appendix [Bl only with the infinitesimals changed to be finite Au(v). G« in Eq. (1 1) of I 
is a 2 V density so that, even though the factor of 1/2 there is neglected, the summation over k of it 
does not equal to Gj, which is IV, in the left column on that same page. [Gj(k) there should have 
been used to denote the IV discrete density instead of the 2V density c$ *(k) in Eq. (11) there. For 
this historical reason, we resist using Gi in this paper.] 

Note adding an arbitrary term, which vanishes as the lattice spacing goes to zero, will not 
change the convergence to the density [of Eq. ([8]), or of Eq. © after transformed to the IV density 
- see Appendix HH with continuum velocity. The additional term should be determined by other 
condition(s). In our case, the "additional" term comes from the two dimensional (covariant) den- 
sity which is determined by the conserved quantities. For our 2V covariance function in continuum 
freedoms of u and v, there is also a second term, the discretization of which exactly corresponds 
to the "additional" term for the discrete density, besides the one with the Dirac delta function. The 
effects of this "additional" term as discussed in Sec. [jV]have been partly confirmed by ongoing 
numerical work. 10 One should also be careful that (g(k,Vi)g*(k,Vi)) = c i;i (k) here for the dis- 
cretized version is not equivalent to Eq. © evaluated at u = v = vt, but to the coarse-graining 
of it [g(k, Vi) is the solution to the discretized system but not <?(k, v) (the solution of the original 
system, evaluated at v = ViJ\. 

We thus have been able to go back and forth smoothly between the distribution over continuous 
u, v (the zero lattice spacing limit) and the discretized version. Eqs. (fT2l) and (fT4l) are the same as 
those obtained earlier in I where the calculations were performed with the discretized system of 
Eqs. (OQ) and © as the starting point. To take the continuum limit, there may be many other ways 
which do not lead to the convergence to our results with continuum degrees of freedom or do not 
have any relevance to the discretization scheme. Some particularly strange limit was indicated in 
Appendix D of I; see Sec. |V]for further remarks. We will not discuss any of them in detail - the 
purpose of this paper is to present the scheme for bridging the continuum-degree density and the 
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discretized one in a consistent way for both 2D and 3D (see below) gyrokinetics. The temperatures 
a(v) and fto are determined by G(v) and E and should not change too much with the discretization 
scheme and resolution, and especially a(vi) and a should converge to a(uj) and a . 



B. Partition function 

The gyrokinetic system was discretized from the beginning in I, however here we will start 
from the discretization of Eqs. dH |5]), as often done for the Feynman path integral. We will 
show how the same results can be obtained and how the continuum limit is recovered. We now 
denote g(lt,Vi) by g kji for simplicity. Using Eq. (fTTI) to discretize Eqs. dUEL we get the Gibbs 
distribution, the same as in I, 

n k 5 k x exp j- ^[Sijai + a 27r/3(k)w i (/c)w i (A;)]^ M ^j| ) 

with the partition function for each k 

z k = n N dst{6ijOii + a 27r/3(k)w i (A;)w j (A;)}" 1 . (15) 

[A factor of 2 N is eliminated by the realizability condition as in the previous subsection.] The term 

$ij&igk,ig)i,j comes from 
6' 

— a(fi)jk,i?k jTrnrrij <= 5(u- v)a(v)g(k,u)g*(k,v)dudv, 
rrij ,J 

which exactly shows how the discretization scheme works. Invoking the matrix determinant 
lemma,— we have 

det{5ijai + ao2nf3(k)wi(k)wj(k)} = [1 + 5ija^ 1 a 2nf3{k)wi(k)wj(k)] det{5i 7 jOZi} 

= [1 + y^^a,~ 1 a 27r/3(k)w i (/;:)w j (fc)]ri i aj. 

i,i 

Substituting this into Eq. (fT5T) . we get the partition function. (It is possible to proceed from here 
to obtain the functional determinant and the functional partition function: Formal divergence may 
be met but shall not cause essential difficulty, as what we really need are the derivatives.) Then, 
we can derive the observables (by introducing the auxiliary function, if necessary). In particular, 



d In z k 1 



ft; 



a 2TTl3(k)w^a 



2„,-l 



l + a 2iTl3(k)J2 l wfa l 



(16) 



which converges to the continuum result through the scheme for discretization/quantization given 
in the previous subsection. This equation is the same as in I. The nice thing is that (g£ j#k,i) is 
exactly conserved by the discretized dynamics as the starting point of I. 

The above brief supplementary material thus consistently shows that the continuum calculation, 
if discretization is postpone to the intermediate stage, yields the same multivariable distribution 
as that from the discretized dynamics. Then, by an alternative method, the partition function, we 
obtain the same final results as I. These two approaches are conceptually different, but consis- 
tent. Especially, the way we obtain the discretized Gibbs distribution highlights the discretization 
scheme through which the results converge to the continuum results. 

IV. DISCUSSIONS 

For the physical discussions here, generally either IV or 2V distribution will work: For the 
case with continuous velocity, obviously the IV density is more convenient; for the quantized case 
with rrii = Av, the corresponding IV distribution is just Eq. (fl4l) multiplied by Av following 
Appendix [B] so that we will just use the 2V density, Eq. (PT4l) . without the necessity of extra work. 

Eqs. ©, © and (flOl ) with negative a seem to argue for the dual-cascade scenario as sug- 
gested by Kraichnan in 2D Navier-Stokes turbulence. 12 The idea is simply that the system persists 
in relaxing to the absolute equilibrium state, with E and G concentrating more at, respectively, 
low- and high-A; modes (some minor oscillations may be introduced by J , but they should not 
change the picture essentially). However, the equilibrium state is continuously destroyed by the 
dissipation (, or other sink, and pumping). When there is a minimal cutoff scale in velocity space 
as introduced by the discretization of velocity shown in Eq. (PT4l) . an extra term emerges. As in 
the discussion of Appendix El but with finite, instead of infinitesimal, Av, the second term in 
Eq. (fT4l) . — ao27r f/, fc ^% _t , gives vanishing contribution to the distribution over v of G when 
Av — > 0. New physics enter when the contribution from this second term is not negligible, which 
happens especially clearly in either of the following two cases where Av is not small enough 
and that the second term can even dominate: The denominator approaches zero from above with 
«o < and ctj > with j > 0, when k — > k m i n ; or, the denominator approaches zero from below 
with olq > and atj < with j > 0, when k — > k max . In general, we have N + 1 conserved 
quantities (totally N lattice points are introduced for the discretization of g on v) defined by the 
summation over k of Eqs. (PT4l) and the discretized version of (flOl) . with their relations constrained 
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by the discretized version of Eq. ©. The spectral relaxation of all these N + 1 conserved quanti- 
ties with the nested constraints would be complicated and it is challenging to specify which case(s) 
is (are) physically relevant to turbulence (note that negative a(vi) = aii/mi is realizable.) We can 
see that the discretized spectral relation between D(k) and Q ^k) could drastically differ from 
that of the continuous v case when a is negative. Especially, when k is small, say, well below 
some injection wavenumber k in , the discretized version of Eq. (flOl) . dominated by the second 
term, may have the same asymptotic spectral behavior as Eq. (fl4l) . In such case the negative-a 
state with most of E and G located at the lowest modes indicates a possibility that all of them are 
transferred to lowest modes simultaneously as one whole conserved quantity. It is also possible 
that finite Larmor effect at intermediate scales may hinder the forward transfer of G so that a large 
amount of G may reside in the low mode regime. Suppose the wavenumber k in where both E and 
G are injected is some intermediate wavenumber, a possible picture would be that both E and G 
are transferred downward to ever larger scales but a forward transfer of G dominates above k in 
to ever smaller scales where dissipation happens. The finite Av introduces a particular truncation 
of velocity scale and can be related to collision/dissipation scale. The results also highlight the 
importance of sufficient numerical resolution in v and the use of appropriate dissipation model or 
collision operator. The former is well known, but to our knowledge, this is the first time the resolu- 
tion effects are calculated in the case of absolute equilibrium. For the latter, a particularly relevant 
and illuminating case from fluid turbulence was shown in Ref. \6. The difference of a factor of Av 
between the IV and 2 V density spectra should also be taken into account in the evaluation of the 
pseudo-dissipation/collision from the thermalization of increasing modes during the transit regime 
of the simulations as done in fluids by Cichowlas et al.— 

With exactly the same ideas and techniques as in the calculations of the 2D case, but extending 
v to v = (v±,vn) (c.f., the end of Appendix I A 2D one can obtain the 3D results which are in 
consistent with Zhu and Hammett, 3 where there was no disagreement between the coauthors in 
taking the continuum limit and discussing the relevant physics. 

V. CONCLUDING REMARKS 

In summary, the exact gyrokinetic statistical absolute equilibria can be obtained analytically 
in two ways: either starting starting from the Galerkin truncated gyrokinetic equations or from 
their discretized versions (with respect to velocity field). That is, the results derived from the 
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original integro-differential equations and those obtained by taking the continuum limit of the 
absolute equilibria of the discretized system are shown to be consistent. The comparison of the 
results brings physical insights such as possible variations of the generic feature of dual cascade 
in magnetized plasmas as a kind of kinetic effect. Most of the physical pictures discussed earlier 
in I with respect to the discretized results, persist in the continuous limit. Therefore, we have not 
repeated them here, but we have only isolated some specific issues. 

The prescription of bridging the absolute equilibria between continuous velocity and quan- 
tized/discretized velocity, 

8 id /mi^5(u - v); ai = a(vi)mu a(vi) a(i>j); <S ->a , 

should also be used for the evaluation of final equilibrium spectra with given initial G(v) and E. 
Specific attention should be paid to the subtle differences between the IV or 2V densities: For 
example, in Appendix D of I, Gi was used for IV densities and subsequently 2V ones in the right 
column of page 10, which may be the reason leading to the final unphysical conclusion, namely, 
in the continuum limit all the energy will be frozen in the lowest mode. [Misconceptions there 
are evidenced by inconsistent conclusions which can be summarized as follows: For small E and 
Av, both E oc k and E(k) = for all k > k min were derived.] The prescription given here for 
regularization also means that more lattice sites (denoted by "i" in Gi there) should appear with 
decreasing lattice spacing, which should have been taken into account in the end of Appendix D 
in I. Gi, evaluated at given velocities (v h Vi), in the main text of I is conserved by the dynamics, 
which however does not mean that it should be fixed when taking the continuum limit: G^s in 
Appendix D of I were fixed as Av — > 0, then either all lattice points were squeezed together (in 
which case % looses its meaning of velocity index) with the total grid number N fixed or an infinite 
number of additional lattice points are inserted, in which case the information of the inserted G 
is undefined. Neither case is physical. The unchanged quantity is the IV density approximated 
by GiAv for Av — > when G is interpolated with given IV density and with given (discretized) 
integral [summation of Gi(Av) 2 ] of this IV density over any velocity interval(s). This is similar 
to many other thermodynamic limit issues where the way of taking the limits is crucial to obtain 
correct or meaningful results. 

This work also serves to outline a procedure for quantifying the errors due to finite truncation 
of velocity scale, although most physically realistic simulations contain collision physics or other 
sources and/or sinks. The methodology may be used to explore other systems with similar features, 
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such as other kinetic models of plasma dynamics, and it can be of interest to researchers concerned 
with developing turbulence model (numerical and analytical) for plasma descriptions based on the 
gyrokinetic ordering. 
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Appendix A: Functional inversion 

Here we sketch the functional inversion problem and related calculations. 

1. A direct approach 

As j 5(u — x)5(x — v)dx = 5(u — v), we can write C(k, u, v) = 2{ v ) +C(k,u,v). Balancing 
the rest terms from Eq. ©, we have 




+ 



a 27r f3 (k)v J (kv)xJ (kx)C (k , u, x)dx. 



a(u) 



1 



/ 



ao27i (3(k)uJo(ku)x Jo(kx)C (k , x, v)dx 



(Al) 



a(u)a(v) a(v) 
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The above formula demonstrates the symmetry of C(k, u, v) and that we can write C(k, u, v) 
c(k, u)c(k, v)c(k), which, combined with Eqs. (I61 IA1I) . gives the equation 

%k)c(k,u)\c{k,v)- VJo W^ m fxMkx)c(k,x)dx] = ^muUku)vHkv) 



a\v) 



a(u)a(v) 



and a similar equation obtained by permuting u and v. We determine (up to a constant which can 

iJo(ku 
a(u) 



be absorbed by c) the form of c as c(k, u) = uJo ^ an d then determine c, and finally C as given 



in Eq. ©. 

2. A general functional inversion formula 

Eq. © can be straight forwardly obtained from the functional inversion formula, 

[A(u,v) + f(u)g(v)] x = [ A x (u,x)[S{x-v) + J f(x)g(y)A x (y,v)dy] I dx 

A T (u, x){5{x - v) - J f(x)g(y)A I (y, v)dy + 
l[j f(x)g(y)A(y,z)dy j f(z)g(y)A I (y,v)dy]dz-...}dx 
A x (u,v) — j A x (u,x)f(x)dx J A x (x, v )g{x)dx\l — JJ f(y)A(x,y)g(x)dxdy + 

+ // f(y)M x ,y)9(x)dxdy // f(y)A(x,y)g(x)dxdy- ...} 



= A x {u v) - t AX((U ' x ^^ dx I Al ( x ' v )9{x)dx (A2) 
1 + 1 1 f(y) AX (x,y)g(x)dxdy 

the discretization of which is the Sherman-Morrison formula for matrix inversion. In the first 
equality of the above formula, we have used [J A(u,w)B(w,v)dw] x = J B x (u, w) A 1 (w , v)dw 
which can be verified directly. The simple formal calculation to obtain Eq. (IA2I) corresponds 
also to the direct extension of Ref. |l4j. To our knowledge, this functional inversion formula 
(where the denominator is of course assumed to be nonzero) in explicit form is new, though it can 
be considered as an extension of the already well known Sherman-Morrion formula for matrix 
inversion. 

u and v in Formula (IA2I) can be vectors in the (3 + 2)-D case. 

Appendix B: Reverting to IV description 

In the main text, we have raised the V-dimension of G in Eqs. (J4l O to 2V. We need to revert 
it to IV. 
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No confusion should be caused as Eq. © involves the distribution on u-v plane while Eq. © 
is the distribution over the coordinate v. Although such reverting to one dimensional distribution 
is intuitively direct, let us consider the following subtleties. The results from Eqs. © and © can 
be understood by invoking the notion of infinitesimals. What we should have is the covariance 
(g(k, u)dug*(k, v )dv) and that we should call (g(k,v)g*(k,v)) the self-correlation density. We 
have 

hm(g{k,u)dug*{k,v)dv) = — — = — 2r2 ' -, (dv) 2 , (Bl) 



«(«) a 2 ( W )[l+2vrao/3(fc)/^S^^ 

where we have applied lim u ^ 5(t> — u)du = l, which can be easily understood with a particular 
nascent delta function rj € (u — v) = e _1 [sgn(e/2 — \u — v\) + l]/2 8{u — v) to regularize the Dirac 
delta and du ^ e, in the way consistent with the prescription for numerical discretization or for 
taking the continuum limit as shown in the main text. Eq. (|Bll) clearly shows that the IV density 
is l/ct(v) which, in nonstandard analysis, should be the standard part of the IV density. [The 
corresponding IV discrete distribution of the 2V discrete distribution Q j can be easily obtained in 
the same way with simply the infinitesimals du(v) changed to finite Au(v) and that an "additional" 
term survives from the second term of Eq. (fT2~l) or CHI)-] The fact is that the second term has 
infinitesimal contribution to the infinitesimal (zigzagged) strip around the line u — v in the u-v 
plane, while the first term with Dirac delta function can be integrated throughout the strip to give 
the finite ID distribution l/a(v): For any reasonably (well) behaved test function T(u, v), 
f f f rv+Au/2 _ n 

/ / C(k,u,v)T(u,v)dudv = lim dv [ ^ h C(k, u, v)]T(u, v)du 

J Jlimu^v Au^Oj J v - Au/2 a(v) 

= [ -^-T(v,v)dv+ lim / {C(k, v, v)T(v, v)Au + 0[(Au) 2 ]}dv = [ -^T(v, v)dv.(B2) 
J oi(v) Au-^oJ J a{v) 

(Interchangability of integration and taking limit is assumed in the above.) 
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